%% Risk Neutral Marke Maker
clc;
clear;

N_M=3;
N_I=10;


lambda_I =1;


T=21;
tau_V=1;



tau_eps=0.5*ones(T,1);



tau_eta=10*ones(T,1);





theta_bar=1;
V_bar=1;   

X_0=0;


theta_I_bar=theta_bar/(N_I+N_M);




o_I_V=zeros(T,1);
for j=1:T
    o_I_V(j)=1/(tau_V+sum(tau_eps(1:j)));
end
K_I_V=o_I_V.*tau_eps;

sigma_I=zeros(T,1);     
sigma_I(1)=1/tau_V+1/tau_eps(1); 
for j=2:T
    sigma_I(j)=o_I_V(j-1)+1/tau_eps(j);  
end

Sigma_I=zeros(2,2,T);
for t=1:T
    Sigma_I(1,1,t)=sigma_I(t);
    Sigma_I(2,2,t)=1/tau_eta(t);
end





mu= lambda_I*o_I_V;


h=lambda_I./tau_eps;

Sigma_U=zeros(T,1);
Sigma_U(1)=1/tau_V+1/tau_eps(1)+h(1)^2/tau_eta(1);

o_U_V=zeros(T,1);
o_U_X=zeros(T,1);
o_U_VX=zeros(T,1);

o_U_V(1)=(tau_eps(1)^-1+h(1)^2*tau_eta(1)^-1)/(tau_V*Sigma_U(1));
o_U_X(1)=1/tau_eta(1)-h(1)^2/(tau_eta(1)^2*Sigma_U(1));
o_U_VX(1)=(1/tau_V)*(h(1)/tau_eta(1))/Sigma_U(1);

K_U_V=zeros(T,1);
K_U_X=zeros(T,1);

K_U_V(1)=tau_V^-1/Sigma_U(1);
K_U_X(1)=-h(1)*tau_eta(1)^-1/Sigma_U(1);

for j=2:T
    Sigma_U(j)=o_U_V(j-1)+1/tau_eps(j)+h(j)^2/tau_eta(j);
    
    o_U_V(j)=(1/tau_eps(j)+h(j)^2/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    o_U_X(j)=o_U_X(j-1)+1/tau_eta(j)-(o_U_VX(j-1)-h(j)/tau_eta(j))^2/Sigma_U(j);
    o_U_VX(j)=o_U_VX(j-1)-(o_U_VX(j-1)-h(j)/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    
    K_U_V(j)=o_U_V(j-1)/Sigma_U(j);
    K_U_X(j)=(o_U_VX(j-1)-h(j)/tau_eta(j))/Sigma_U(j);
end

H_I=zeros(4,4,T);
for t=2:T
    H_I(1,1,t)=1;
    H_I(2,2,t)=1;
    H_I(3,3,t)=1-K_U_X(t)*mu(t-1);
    H_I(3,2,t)=K_U_X(t)*mu(t-1);
    H_I(4,4,t)=1;
end

F_I=zeros(4,2,T);
for t=1:T
    F_I(1,1,t)=K_I_V(t);
    F_I(2,2,t)=1;
    F_I(3,1,t)=K_U_X(t);
    F_I(3,2,t)=-K_U_X(t)*h(t);
end

F_M=zeros(3,1,T);
for t=1:T
    F_M(1,1,t)=K_U_V(t);
    F_M(2,1,t)=K_U_X(t);
end








% the values of the parameters for the last period

gamma_I=zeros(T,1);
gamma_I(T)=lambda_I*o_I_V(T);


delta_I=zeros(T,1);
delta_I(T)=1;

delta_M=zeros(T,1);
delta_M(T)=0;

g_I=zeros(T,1);
g_I(T)=N_M*lambda_I*o_I_V(T)/(N_M+1);

g_M=zeros(T,1);
g_M(T)=-N_I/(N_M+1);

f_I=zeros(T,1);
f_I(T)=lambda_I*o_I_V(T)/(N_M+1);

f_M=zeros(T,1);
f_M(T)=N_I/(N_M+1);

rho_I=zeros(T,1);
rho_I(T)=1;

rho_M=zeros(T,1);
rho_M(T)=0;

m_I=zeros(T,1);
m_I(T)=lambda_I*o_I_V(T)/(N_M+1)^2;

m_M=zeros(T,1);
m_M(T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;




C_I=zeros(4,1,T);
C_I(3,1,T)=-lambda_I*o_I_V(T)*N_M/(N_M+1)^2;


C_M=zeros(3,1,T);
C_M(2,1,T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;

M_I=zeros(4,4,T);
M_I(1,2,T)=1;
M_I(2,1,T)=1;
M_I(2,2,T)=-lambda_I*o_I_V(T);
M_I(2,4,T)=-1;
M_I(3,3,T)=lambda_I*o_I_V(T)*(N_M/(N_M+1))^2;
M_I(4,2,T)=-1;

M_M=zeros(3,3,T);
M_M(2,2,T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;



Xi_I=zeros(2,2,T);
Xi_I(:,:,T)=eye(2)/( eye(2)/Sigma_I(:,:,T)+ lambda_I*transpose(F_I(:,:,T))*M_I(:,:,T)*F_I(:,:,T));



H_I_P=zeros(4,1,T);
H_I_P(1,1,T)=delta_I(T);
H_I_P(2,1,T)=-mu(T)*delta_I(T);
H_I_P(3,1,T)=g_I(T);
H_I_P(4,1,T)=1-delta_I(T);

H_M_alpha=zeros(3,1,T);
H_M_alpha(1,1,T)= -delta_M(T);
H_M_alpha(2,1,T)= mu(T)*delta_M(T)+g_M(T);
H_M_alpha(3,1,T)= delta_M(T);

H_I_alpha=zeros(4,1,T);
H_I_alpha(1,1,T)= -delta_M(T);
H_I_alpha(2,1,T)= delta_M(T)*mu(T);
H_I_alpha(3,1,T)= g_M(T);
H_I_alpha(4,1,T)= delta_M(T);

H_I_theta=zeros(4,1,T);
H_I_theta(:,:,T)=H_I_alpha(:,:,T)*N_M/N_I;

f_I_theta=zeros(T,1);
f_I_theta(T)=1-f_M(T)*N_M/N_I;




H_I_D=zeros(4,1,T);







test1=zeros(T,1);
test1(T)=mu(T);

test2=zeros(T,1);
test2(T)=1;

test3=zeros(T,1);


SOD_I=zeros(T-1,1);
SOD_M=zeros(T-1,1);


% compute the values of the parameters recursively

for t=T:-1:2
    
    
    H_I_D(:,:,t)=transpose( transpose(H_I_P(:,:,t))*H_I(:,:,t)-...
        lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)...
        *transpose(F_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t) );
    
    
    test1(t-1)=-H_I_D(2,1,t)/H_I_D(1,1,t);   % =mu(t-1)
    test2(t-1)=H_I_D(1,1,t)+H_I_D(4,1,t);  % =1
    
       
    gamma_I(t-1)=f_I(t)+lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)...
        *transpose(F_I(:,:,t))*(H_I_P(:,:,t)+C_I(:,:,t));
    
    SOD_I(t-1)=lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*H_I_P(:,:,t);
    SOD_M(t-1)=gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2;
    
    
    delta_M(t-1)=( 1-H_I_D(1,1,t)-C_M(1,1,t)/N_I )/ (gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
    
    g_M(t-1)=( H_I_D(3,1,t)+C_M(2,1,t)/N_I+ mu(t-1)*C_M(1,1,t)/N_I-mu(t-1) )...
        /(gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
    
    f_M(t-1)=( gamma_I(t-1)-m_M(t)/N_I )/(gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
   
    f_I_theta(t-1)=1-f_M(t-1)*N_M/N_I;
    
    H_I_alpha(1,1,t-1)= -delta_M(t-1);
    H_I_alpha(2,1,t-1)= delta_M(t-1)*mu(t-1);
    H_I_alpha(3,1,t-1)= g_M(t-1);
    H_I_alpha(4,1,t-1)= delta_M(t-1);
    
    H_M_alpha(1,1,t-1)=-delta_M(t-1);
    H_M_alpha(2,1,t-1)=delta_M(t-1)*mu(t-1)+g_M(t-1);
    H_M_alpha(3,1,t-1)=delta_M(t-1);

    H_I_theta(:,:,t-1)=H_I_alpha(:,:,t-1)*N_M/N_I;
    
    delta_I(t-1)=H_I_D(1,1,t)+gamma_I(t-1)*delta_M(t-1)*N_M/N_I;
    
    g_I(t-1)=H_I_D(3,1,t)-gamma_I(t-1)*g_M(t-1)*N_M/N_I;
    
    f_I(t-1)=gamma_I(t-1)-gamma_I(t-1)*f_M(t-1)*N_M/N_I;
    
    H_I_P(1,1,t-1)=delta_I(t-1);
    H_I_P(2,1,t-1)=-mu(t-1)*delta_I(t-1);
    H_I_P(3,1,t-1)=g_I(t-1);
    H_I_P(4,1,t-1)=1-delta_I(t-1);
    
    
    
    
    
    m_I(t-1)= (f_I_theta(t-1))^2* (m_I(t)+lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- lambda_I*transpose(C_I(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t) );
      
    
    m_M(t-1)=m_M(t)*(1-2*f_M(t-1)*(N_M-1)/N_I) +...
        f_M(t-1)^2* (2*gamma_I(t-1)/N_I+ m_M(t)* (N_M-2)*N_M/(N_I)^2);
    
   
    
    C_I(:,:,t-1)=f_I_theta(t-1)*( transpose(H_I(:,:,t))*C_I(:,:,t)+ m_I(t)*H_I_theta(:,:,t-1)+...
        lambda_I* H_I_theta(:,:,t-1)*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- lambda_I*f_I_theta(t-1)*transpose( C_I(:,:,t)*...
        transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) )*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t) );
    
    C_M(:,:,t-1)= (1-f_M(t-1)*(N_M-1)/N_I)*C_M(:,:,t)+...
        H_M_alpha(:,:,t-1)*...
        (m_M(t)*(N_M-1)/N_I- m_M(t)*f_M(t-1)*(N_M-2)*N_M/(N_I)^2-2*f_M(t-1)*gamma_I(t-1)/N_I );
    
    test3(t-1)=C_M(1,1,t-1)+C_M(3,1,t-1); %=0
    
    
    
    M_I(:,:,t-1)=transpose(H_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t)+m_I(t)*H_I_theta(:,:,t-1)*...
        transpose(H_I_theta(:,:,t-1))+ H_I_theta(:,:,t-1)*transpose(C_I(:,:,t))*H_I(:,:,t)+...
        transpose(H_I_theta(:,:,t-1)*transpose(C_I(:,:,t))*H_I(:,:,t))+...
        lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*...
        H_I_P(:,:,t)*H_I_theta(:,:,t-1)*transpose(H_I_theta(:,:,t-1))-...
        lambda_I*transpose( C_I(:,:,t)*transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) )*...
        F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*( C_I(:,:,t)*...
        transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) );
    
    M_M(:,:,t-1)=M_M(:,:,t)+ (2*gamma_I(t-1)/N_I+ m_M(t)*(N_M-2)*N_M/(N_I)^2)*H_M_alpha(:,:,t-1)*...
        transpose(H_M_alpha(:,:,t-1)) + H_M_alpha(:,:,t-1)*transpose(C_M(:,:,t))*(N_M-1)/N_I+...
        C_M(:,:,t)*transpose(H_M_alpha(:,:,t-1))*(N_M-1)/N_I;
    
    
    
    Xi_I(:,:,t-1)= eye(2)/(eye(2)/Sigma_I(:,:,t-1)+...
        lambda_I*transpose(F_I(:,:,t-1))*M_I(:,:,t-1)*F_I(:,:,t-1));

    
    rho_I(t-1)=rho_I(t)*sqrt(det(Xi_I(:,:,t))/det(Sigma_I(:,:,t)) );
    
    rho_M(t-1)=rho_M(t)+ transpose(F_M(:,:,t))*M_M(:,:,t)*F_M(:,:,t)*Sigma_U(t)/2;   

end

rho_M_0=rho_M(1)+ transpose(F_M(:,:,1))*M_M(:,:,1)*F_M(:,:,1)*Sigma_U(t)/2; 

SOC_I=find(SOD_I<0)
SOC_M=find(SOD_M<0)








mu_alpha=zeros(T,1);
mu_alpha(1)=-f_M(1)*theta_I_bar;


mu_D=zeros(T,1);
mu_D(1)=-f_I(1)*theta_I_bar;




for t=2:T
    mu_alpha(t)=-f_M(t)*theta_I_bar*prod(f_I_theta(1:t-1));
    mu_D(t)=-f_I(t)*theta_I_bar*prod(f_I_theta(1:t-1));
end


    











% Expected trading volume and bid-ask spread

tau=(1./tau_eps+h.^2./tau_eta).^-1;



d_V=zeros(T,1);
d_eps=zeros(T,T);
d_eta=zeros(T,T);

k_V=zeros(T,1);
k_eps=zeros(T,T);
k_eta=zeros(T,T);



for t=1:T
    
    d_V(t)= (g_M(t)/mu(t))*o_U_V(t)*sum(tau(1:t)) -(delta_M(t)+g_M(t)/mu(t))*o_I_V(t)*...
        sum(tau_eps(1:t));
    k_V(t)=-(1-g_I(t)/mu(t))*o_U_V(t)*sum(tau(1:t))+ (delta_I(t)-g_I(t)/mu(t))*o_I_V(t)*...
        sum(tau_eps(1:t));
    
    for n=1:t
        d_eps(t,n)= (g_M(t)/mu(t))*o_U_V(t)*tau(n)- (delta_M(t)+g_M(t)/mu(t))*o_I_V(t)*tau_eps(n);
        k_eps(t,n)=-(1-g_I(t)/mu(t))*o_U_V(t)*tau(n)+ (delta_I(t)-g_I(t)/mu(t))*o_I_V(t)*tau_eps(n);
       
        
        d_eta(t,n)= delta_M(t)*mu(t)+ g_M(t)- (g_M(t)/mu(t))*o_U_V(t)*tau(n)*h(n);
        k_eta(t,n)=-(delta_I(t)*mu(t)-g_I(t))+(1-g_I(t)/mu(t))*o_U_V(t)*tau(n)*h(n);
    end
end




D_V=zeros(T,1);
D_V(1)=d_V(1);
D_V(2)=d_V(2)-(N_M/N_I)*f_M(2)*d_V(1);

K_V=zeros(T,1);
K_V(1)=k_V(1);
K_V(2)=k_V(2)-(N_M/N_I)*f_I(2)*d_V(1);

for t=3:T
    
    c_V=zeros(t-1);
    for i=2:t-1
        c_V(i)=prod(f_I_theta(t+1-i:t-1))*d_V(t-i);
    end 
    
    D_V(t)=d_V(t)-(N_M/N_I)*f_M(t)* (d_V(t-1)+ sum(c_V(2:t-1)) );
    K_V(t)=k_V(t)-(N_M/N_I)*f_I(t)* (d_V(t-1)+ sum(c_V(2:t-1)) );
     
end

D_eps=zeros(T,T);
D_eps(1,1)=d_eps(1,1);


K_eps=zeros(T,T);
K_eps(1,1)=k_eps(1,1);


D_eta=zeros(T,T);
D_eta(1,1)=d_eta(1,1);


K_eta=zeros(T,T);
K_eta(1,1)=k_eta(1,1);



for t=2:T
    D_eps(t,t)=d_eps(t,t);
    D_eps(t,t-1)=d_eps(t,t-1)-(N_M/N_I)*f_M(t)*d_eps(t-1,t-1);
    
    K_eps(t,t)=k_eps(t,t);
    K_eps(t,t-1)=k_eps(t,t-1)-(N_M/N_I)*f_I(t)*d_eps(t-1,t-1);
    
    D_eta(t,t)=d_eta(t,t);
    D_eta(t,t-1)=d_eta(t,t-1)-(N_M/N_I)*f_M(t)*d_eta(t-1,t-1);
    
    K_eta(t,t)=k_eta(t,t);
    K_eta(t,t-1)=k_eta(t,t-1)-(N_M/N_I)*f_I(t)*d_eta(t-1,t-1);
    
    
    for n=1:t-2
        c_eps=zeros(t-n);
        c_eta=zeros(t-n);
        for i=2:t-n
            c_eps(i)=prod(f_I_theta(t+1-i:t-1))*d_eps(t-i,n);
            c_eta(i)=prod(f_I_theta(t+1-i:t-1))*d_eta(t-i,n);
        end
        D_eps(t,n)=d_eps(t,n)-(N_M/N_I)*f_M(t)*d_eps(t-1,n)-(N_M/N_I)*f_M(t)*sum(c_eps(2:t-n));
        K_eps(t,n)=k_eps(t,n)-(N_M/N_I)*f_I(t)*d_eps(t-1,n)-(N_M/N_I)*f_I(t)*sum(c_eps(2:t-n)); 
        
        D_eta(t,n)=d_eta(t,n)-(N_M/N_I)*f_M(t)*d_eta(t-1,n)-(N_M/N_I)*f_M(t)*sum(c_eta(2:t-n));
        K_eta(t,n)=k_eta(t,n)-(N_M/N_I)*f_I(t)*d_eta(t-1,n)-(N_M/N_I)*f_I(t)*sum(c_eta(2:t-n));
    end
end


            

sigma_alpha_sq=zeros(T,1);
for t=1:T
    sigma_alpha_sq(t)= D_V(t)^2/tau_V + sum(D_eps(t,:).^2./transpose(tau_eps))+...
        sum(D_eta(t,:).^2./transpose(tau_eta));
end
sigma_alpha=sqrt(sigma_alpha_sq);

sigma_D_sq=zeros(T,1);
for t=1:T
    sigma_D_sq(t)= K_V(t)^2/tau_V + sum(K_eps(t,:).^2./transpose(tau_eps))+...
        sum(K_eta(t,:).^2./transpose(tau_eta));
end
sigma_D=sqrt(sigma_D_sq);





Vol= N_M*(  sqrt(2/pi)*sigma_alpha.*exp( -0.5*(mu_alpha./sigma_alpha).^2 ) +...
    mu_alpha.*( 1-2* normcdf(-mu_alpha./sigma_alpha) ) ) ;




Spread=  sqrt(2/pi)*sigma_D.*exp( -0.5*(mu_D./sigma_D).^2 ) +...
    mu_D.*( 1-2* normcdf(-mu_D./sigma_D) )  ;










figure(1)
plot(0:1/(T-1):1, Vol);
xlabel('$t/T$','interpreter','latex')
ylabel('Expected Trading Volume','interpreter','latex')



figure(2)
plot(0:1/(T-1):1, Spread, 'r');
xlabel('$t/T$','interpreter','latex')
ylabel('Expected Bid-Ask Spread','interpreter','latex')








    



